Chapter 18:- NumPy, SciPy From book Python Programming (Problem solving, Packages and Libraries) published by McGraw Hill Education (India) Private limited. By:

  • Anurag Gupta
  • G. P. Biswas

Note the following:-

  1. This html document is meant as an accompaniment to Chapter 18 NumPy, SciPy .
  2. The document contains scripts executed on IDLE as well as on Jupyter notebook.
  3. The scripts executed on Jupyter can be directly copied and run into a Jupyter notebook or some other IDE (Like Pycharm or Eclipse with PyDev or Visual studio).
  4. However the scripts on IDLE also contain the >>> symbol and therefore cannot be directly executed. If you want to execute them on IDLE or Jupyter, you need to manually remove the >>> symbol.
  5. Wherever needed some background material from the book is also included to help you better understand the scripts
  6. The topic numbers given on each paragraph match the topic numbers of the book, so you can easily identify the topics and corresponding scripts.
  7. In some of the scripts, the file paths give are that of the author's computer. You need to replace them with file paths of your own computer.
  8. The topic on "broadcasting of arrays" contains some additional material which is not given in the book.
  9. At some places, to improve readability, page numbers of the book are indicated in green font like:- See Page 181 of the book
  10. This document was first created as a Jupyter Notebook as combination of Markdown and code cells (extension .ipynb) and then downloaded as html. If someone wants to "modify" or "extend' this document, you may ask for the original .ipynb file by sending me an e-mail at:- 999.anuraggupta@gmail.com

18.2. Basics of NumPy and SciPy
18.2.1. N-dimensional array in NumPy

The main data structure defined and used in NumPy module is called ndarray. The term ndarray stands for N-dimensional array. An object of type ndarray represents a multi-dimensional array of data.

The use of NumPy will become clear from the following script:
This script is available on page 448 of the book

In [1]:
import numpy as np
A = [ n * n for n in range(20)]
B = np.array(A)
print(B)
C = B * 2 #You cant do this in Python but you can in numpy
print(C)
import matplotlib.pyplot as plt
plt.plot(C)
plt.show()

[  0   1   4   9  16  25  36  49  64  81 100 121 144 169 196 225 256 289
 324 361]
[  0   2   8  18  32  50  72  98 128 162 200 242 288 338 392 450 512 578
 648 722]

However, the data in an array of type ndarray should be of uniform type. The ndarray object determines the type of data at the time of construction of the ndarray object. However, this does not mean that a NumPy array cannot store data of heterogeneous types.

NumPy wraps each heterogeneous object item in the array into a common object format. The array() constructor of a NumPy array has an attribute dtype which must be set to object in case you use different types of objects in a NumPy array. This will be clear from the following code:
This script is available on page 449 of the book

In [2]:
import numpy as np
A = np.array([1, 'cat', True, [1,2]], dtype = object) # Note use of argument dtype
print('Array is ->', A)
print('Array is of type->', A.dtype)

Array is -> [1 'cat' True list([1, 2])]
Array is of type-> object

18.2.2. Some NumPy methods and properties

The NumPy package has many methods and properties to create and manipulate arrays. Some are described in the following text:

NumPy array creation:
This topic is given on page 449 of the book
Arrays can be created with NumPy in a number of ways. Also note that the NumPy package is included in the SciPy package. So all the ways of creating arrays in NumPy can also be done in SciPy. Some of the common methods available in NumPy (and therefore in SciPy) are:

  1. numpy.array(alist), where alist is a python list: It is a function to create a NumPy array which is of type ndarray. Note that in pure Python there are only lists and there is no data type called ndarray. But NumPy has a large number of its own data types. Note that numpy.array() is a method and the object created by this method is ndarray
  2. numpy.zeros(N, M, dtype = float): This constructs an n-dimensional array of the specified shape, i.e., N x M, filled with zeros of the specified dtype.
  3. numpy.ones(N, M, dtype=float): This constructs an 2-dimensional array of the specified shape, i.e., N x M, filled with ones of the specified dtype.
  4. numpy.eye(N, M = None, dtype=float): This returns a two-dimensional (2D) array which has number 1 on the diagonal elements and number 0 everywhere else. By default the data type, i.e., the 0s and the 1s are floats. Note if value of M is not specified, it will be same as N and the resulting matrix will be a square matrix. But if M is different from N then it will not be a square matrix.
  5. numpy.transpose(a, axes=None): This method permutes the dimensions of the given array, i.e., “a”. Note that the parameter a is an array, whose transpose is to be created. You have the option of not specifying any value for axes. If no value is given, then the default value for axes is None and the dimensions get reversed. The other options for axes are not discussed here. You can always look up the online documentation for details. However, note that the axes are passed as a tuple of integers starting from 0 since the 1st axes will have a dimension number of 0. In its simplest form, if you create the transpose of a 2D matrix say of type 2 x 5, then its transpose will be a matrix of type 5 x 2.
  6. numpy.arange(start, stop, increment): NumPy method arrange() is like the Python range() function but with one important difference. In the arrange() method, the values for start, stop and step can be real values (not just integers like in ordinary Python), i.e., can be decimals and fractions also.
  7. numpy.random.random((N, M)): This will create a NumPy array N x M filled with random numbers in range [0,1). Note 1 is not generated but 0 can be generated.
  8. numpy.random.randint(X, Y, (N, M)): Here X and Y must be integers and so random integers will be generated between [X, Y), i.e., from X to Y but not including Y. N  M are the dimensions of the array created.

Examples of use of all the above methods are given in the following code:
This script is available on page 450 of the book

In [3]:
import numpy as np
A = np.array([1, 2, 3, 4]) # Create numpy array from python list
print('A ->', A)
B = np.zeros((2, 2)) # Data type defaults to float
print('B->', B)
C = np.zeros((2,2), dtype = int) # Data type is int, so no decimal after 0
print('C->', C)
D = np.ones((2, 2), dtype = complex) # dtype is complex
print('D->', D)
E = np.eye(4, dtype = int)
print('E->',E)
F = np.array([[1,2], [3,4], [5,6], [7,8]]) # F is 4 x 2 array
print('F->', F)
G = np.transpose(F) # G is 2 x 4 array
print('G->', G)
H = np.arange(2.4, 3.4, 0.2) # start, stop and step are decimals.
print('H->', H)
I = np.random.random((2, 2)) # Create 2 x 2 matrix of randoms in range [0,1)
print('I->', I)
J = np.random.randint(2, 4, (2, 3)) # Will generate random number 2,3 but not 4
print('J->', J)

A -> [1 2 3 4]
B-> [[0. 0.]
 [0. 0.]]
C-> [[0 0]
 [0 0]]
D-> [[1.+0.j 1.+0.j]
 [1.+0.j 1.+0.j]]
E-> [[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]
F-> [[1 2]
 [3 4]
 [5 6]
 [7 8]]
G-> [[1 3 5 7]
 [2 4 6 8]]
H-> [2.4 2.6 2.8 3.  3.2]
I-> [[0.34438222 0.57284245]
 [0.20326575 0.81209942]]
J-> [[3 3 3]
 [2 2 3]]

Dimensions of an array

The concept of dimensions of an array in NumPy will become clear from the following example. As an example, take a $2 \times 3 \times 4$ array. It will have a total of 24 items in it. Fill it up with integers from 0 to 23. The code is as follows:
This script is available on page 451 of the book

In [4]:
import numpy as np
A = np.arange(24).reshape(2,3,4)
print(A)

[[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]]

You may modify the above code to make the concept of dimensions clearer. Consider the following code:
This script is available on page 452 of the book

In [5]:
import numpy as np
A = np.arange(24).reshape(2,3,4)
print('(Dim 0-pos 0)->')
print(A[0, : , :])
print('(Dim0, pos0), (Dim1, pos1)->', A[0, 1, :])
print('(Dim0, pos0), (Dim1, pos1), (Dim2 pos2)->', A[0, 1, 2])

(Dim 0-pos 0)->
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
(Dim0, pos0), (Dim1, pos1)-> [4 5 6 7]
(Dim0, pos0), (Dim1, pos1), (Dim2 pos2)-> 6
Note: In NumPy, axes are defined for those arrays which have more than one dimension. For example a 2D array has two axis. The first axis of a 2D array represents the rows and the second dimension in the array represents the columns in the array. In Mathematics/Physics, you may define dimensions as “The minimum number of coordinates needed to specify a point”. But in NumPy, dimensions is same as axis.

This will become clear from the following example:
This script is available on page 453 of the book

In [6]:
import numpy as np
A = np.array([[10, 20, 30], [40, 50, 60]])# 2 x 3 array
print('shape of A ->', A.shape)
B = A.reshape(3,2) # B will be a 3 x 2 array
print('shape of B->', B.shape)
print('First row->', A[0])
print('Third column->', A[:,2])

shape of A -> (2, 3)
shape of B-> (3, 2)
First row-> [10 20 30]
Third column-> [30 60]

numpy.reshape:

NumPy arrays have a reshape method. This method allows the number of dimensions, and the size of each dimension to be changed. However, the total number of array elements must remain the same. This will be clear from the following examples:
This script is available on page 453 of the book

In [7]:
import numpy as np
A = np.array([[10, 20], [30, 40]])
print('Array A->',A)
print('A as 4 x 1 matrix->', A.reshape(4,1))
print('A as 1 x 4 matrix->', A.reshape(1,4))

Array A-> [[10 20]
 [30 40]]
A as 4 x 1 matrix-> [[10]
 [20]
 [30]
 [40]]
A as 1 x 4 matrix-> [[10 20 30 40]]

Using the repeat method of the NumPy array, you can create an array with repeated elements as shown in the following code:
This script is available on page 454 of the book

In [8]:
import numpy as np
A = np.array([[10, 20], [30, 40]])
B = np.repeat(A, 3)
print(B)

[10 10 10 20 20 20 30 30 30 40 40 40]

Note that the original array A has been converted to a 1D array. To retain the original dimensions, you need to write the script as follows:
This script is available on page 454 of the book

In [9]:
import numpy as np
A = np.array([[10, 20], [30, 40]])
B = np.repeat(A, 3, axis = 0)
print(B)
C = np.repeat(A, 3, axis = 1)
print(C)

[[10 20]
 [10 20]
 [10 20]
 [30 40]
 [30 40]
 [30 40]]
[[10 10 10 20 20 20]
 [30 30 30 40 40 40]]

Element-wise operations on arrays:

In the NumPy package, a number of methods are provided for doing operations on the individual elements of the array. For example, numpy.sum(a) will add all the elements of an array and give the result.
This script is available on page 454 of the book

In [10]:
import numpy as np
A = np.arange(12).reshape(3, 4) # Create 3 x 4 array
print("A=",A)
B = A.sum(axis = 0) # Add along axis 0 
print("B=",B)
C = A.sum(axis = 1) # Add along axis 1
print("C=", C)

A= [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
B= [12 15 18 21]
C= [ 6 22 38]

18.2.4. Broadcasting in NumPy array operations

In NumPy array operations, the term broadcasting comes into play when you do arithmetic operations on arrays with different shapes. In normal Python arrays like lists, mathematical operations are done on arrays of same shape.

However, in NumPy, under certain circumstances, it is possible to do some mathematical operations on arrays of different lengths. The way this is done is to broadcast the smaller array into an array of larger dimension. However considering broadcasting, you need to understand that Python arrays and NumPy arrays behave differently. The following example shows the difference in behavior of Python array (i.e. list) and NumPy array:

In [11]:
import numpy as np
# In Python
myL = [1, 2, 3, 4]
expand = 2
expandedL = myL * 2
print('in python->', expandedL)
# In numpy
a = np.array(myL)
expanded_a = a * expand
print('in numpy->', expanded_a)

in python-> [1, 2, 3, 4, 1, 2, 3, 4]
in numpy-> [2 4 6 8]

Broadcasting of a “scalar”
This topic is not covered in the book
The simplest broadcasting is when an array and a “scalar” are combined in an operation. A “scalar” here means simply a number. The scalar is “converted” into an array of same dimension as the array on which it is being operated and each element in this “converted array” is a copy of the scalar. Before considering broadcasting, of arrays of different shapes, you need to understand that there are 2 aspects of an ndarray:

  1. The number of dimensions or the number of axes of the array and
  2. The size of each dimension. For example if you have an array say 2 x 3 x 4, then it has three dimensions or axes. Further the sizes of the dimensions are 2, 3 and 4.

The rules for broadcasting are:

  • The two arrays have exactly the same shape. (Here of course no broadcasting is required)
  • The two arrays have the same number of dimensions but the sizes of some or all of the dimensions are different. For example you have 2 arrays of dimensions 2 x 3 x 4 and 2 x 3 x 1. These 2 arrays both have 3 dimensions out of which the first 2 dimensions are of same size but the third dimension is of different size. The rule for broadcasting here is that the “uncommon” dimensions should be 1. So for example if you have two 3-D arrays say A → 2 x 3 x 4 and B → 2 x 3 x 1, then array B can be broadcast into 2 x 3 x 4 array because its uncommon dimension has a size of 1.
  • The 2 arrays have different dimensions. It is possible to broadcast an array of smaller dimension into an array of bigger dimension by “Lining up the sizes of the trailing axes of these arrays according to the broadcast rules”. This appears a bit confusing but can be easily understood with examples. Suppose you have array A → 2 x 3 x 4 and array B → 3 x 4. Here the “trailing dimensions” of A and B are matched and both are 3 x 4. So array B is broadcast into an array of 2 x 3 x 4, so as to match the 2 arrays. In this example the “trailing dimensions” matched. However they need not match. Trailing dimensions can be broadcast even if one of the dimensions is 1. So if you had A → 2 x 3 x 4 and B → 1 x 4, then also B can be broadcast into 2 x 3 x 4. This can be done because the “non-matching” dimension of B is 1. Note that the non-matching dimension can be in either of the 2 arrays. So if you have A → 2 x 1 x 4 and B → 3 x 4, then A will be broadcast to A → 2 x 3 x 4 and B will be broadcast to B → 2 x 3 x 4. Here The size of one of the dimensions of A has been increased from 1 to 3, while the number of dimensions of B have been increased from two to three.

Following example will clarify the concepts:

Case1: Multiplying an array with a scalar
This script is not given in the book

In [12]:
import numpy as np
A = np.arange(12).reshape(3, 4) # Create 3 x 4 array
print(A)
scalar = 5
print(A * scalar)

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[[ 0  5 10 15]
 [20 25 30 35]
 [40 45 50 55]]

Note that each element of array A gets multiplierd by the scalar. So here the scalar was broadcast into an array of 3 x 4.

Case2: The 2 arrays are of different dimensions but the “size” of trailing dimensions are same
This script is not given in the book

In [13]:
import numpy as np
# Trailing dimensions of A and B are same i.e 3 x 4
A = np.arange(24).reshape(2, 3, 4) # Create 2 x 3 x 4 array
B = np.arange(12).reshape(3, 4)  # Create a 3 x 4 array

print('A->', A)
print('B->', B)
print('A * B->', A * B)

A-> [[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]]
B-> [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
A * B-> [[[  0   1   4   9]
  [ 16  25  36  49]
  [ 64  81 100 121]]

 [[  0  13  28  45]
  [ 64  85 108 133]
  [160 189 220 253]]]

Case3: The two arrays are of different dimensions and also of different sizes. But the dimension which varies on size has a size of 1.
This script is not given in the book

In [14]:
import numpy as np
A = np.arange(8).reshape(2, 1, 4) # Create 2 x 1 x 4 array
B = np.arange(12).reshape(3, 4)  # Create a 3 x 4 array
# A will be broadcast from (2 x 1 x 4) to (2 x 3 x 4)
# B will be broadcast from (3 x 4) to (2 x 3 x 4)
print('A->', A)
print('B->', B)
print('A * B->', A * B)

A-> [[[0 1 2 3]]

 [[4 5 6 7]]]
B-> [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
A * B-> [[[ 0  1  4  9]
  [ 0  5 12 21]
  [ 0  9 20 33]]

 [[ 0  5 12 21]
  [16 25 36 49]
  [32 45 60 77]]]

18.2.5. Array indexing in NumPy

The use of square brackets ([]) to access items of an array is called array indexing.
For 1D array, indexing in NumPy is exactly like Python list, but for multi-dimensional arrays, there can be a slight difference. This will become from the following example:
This script is available on page 456 of the book

In [15]:
import numpy as np
a = np.arange(12)
a.shape = (4,3)
print(a)
# python way of calling can be used
print(a[-1][-1])
# but in numpy following method can also be used
print(a[-1, -1])  # no need for 2 sets of square brackets

import numpy as np
a = np.arange(12)
a.shape = (4,3)
print(a)
# python way of calling can be used
print(a[-1][-1])
# but in numpy following method can also be used
print(a[-1, -1])  # no need for 2 sets of square brackets

[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
11
11
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
11
11

In Python one can stride over a list using a format similar to my_list[start:end:step]. In NumPy you can do this over a multi-dimensional array shown as follows:
This script is available on page 456 of the book

In [16]:
import numpy as np
a = np.arange(42)
a.shape = (6, 7)
print(a)
b = a[0:6:2, 0:7:3]
print(b)
c = a[::, 0:7:2]
print(c)

[[ 0  1  2  3  4  5  6]
 [ 7  8  9 10 11 12 13]
 [14 15 16 17 18 19 20]
 [21 22 23 24 25 26 27]
 [28 29 30 31 32 33 34]
 [35 36 37 38 39 40 41]]
[[ 0  3  6]
 [14 17 20]
 [28 31 34]]
[[ 0  2  4  6]
 [ 7  9 11 13]
 [14 16 18 20]
 [21 23 25 27]
 [28 30 32 34]
 [35 37 39 41]]

18.2.6. Infinity, negative infinity, zero, negative zero and some other constants in NumPy
In many scientific calculations, you need to represent positive and negative infinity and also positive and negative zero. In addition, certain calculations may not yield a number (for example, log of a negative number is not a number (NaN)). In NumPy NaN is not infinity. Table 18.1 (In the book) gives some common constants in NumPy.
Following script shows usage of these constants:
This script is available on page 457 of the book

In [17]:
import numpy as np
print('Negative infinity->', np.NINF)
print('Positive infinity->', np.PINF)
print('Negative Zero->', np.NZERO)
print('Positive Zero->', np.PZERO)
print('Eulers constant->', np.e)
print('Pi->', np.pi)
print('euler_gamma->', np.euler_gamma)
print('log 0 is infinity->', np.log(0) )
print('log -1 is nan->', np.log(-1) )

Negative infinity-> -inf
Positive infinity-> inf
Negative Zero-> -0.0
Positive Zero-> 0.0
Eulers constant-> 2.718281828459045
Pi-> 3.141592653589793
euler_gamma-> 0.5772156649015329
log 0 is infinity-> -inf
log -1 is nan-> nan

C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:9: RuntimeWarning: divide by zero encountered in log
  if __name__ == '__main__':
C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in log
  # Remove the CWD from sys.path while we load stuff.

18.2.7. np.linspace
NumPy has a method linspace. This method is very useful in getting evenly spaced numbers over an interval. Meaning of evenly spaced numbers:

  • Suppose you have an interval say [2, 4].
  • Suppose you want to divide this interval into 6 intervals then your list would be: [2. 2.4 2.8 3.2 3.6 4. ]. So this is the case where you divided an interval with start = 2, stop = 4 and interval = 6. Note here you also included the endpoint to the right, i.e., 4; so you did endpoint = True.
  • Now normally in many Python functions such as range(), the end point to the right is excluded. (For example, range(10) produces numbers from 0 to 9.) You could do this with the linspace method also by having endpoint = False.
    The signature of linspace method is:
    numpy.linspace(start, stop, interval, endpoint=True)
    This will become clear from the following example:
    This script is available on page 458 of the book
In [18]:
import numpy as np
x1 = np.linspace(2, 4, 6, endpoint = True)
print(x1)
x2 = np.linspace(2, 4, 5, endpoint = False) # Have to reduce endpoint by 1 to get similar intervals
print(x2)

[2.  2.4 2.8 3.2 3.6 4. ]
[2.  2.4 2.8 3.2 3.6]

18.2.8. Understanding np.meshgrid()
To be able to do 3D plotting in Matplotlib (which is explained in later chapter 21), you first need to understand the concept of np.meshgrid() method. Note np is the conventional import for NumPy.
To plot a function in 3D say x, y and z, you could have z = f(x,y) where f(x,y) is some function of x and y. Now both x and y will have some domain and z will have a corresponding range. Suppose the domain of x was integers from 1 to 4. So $x\ ϵ \ (1, 2, 3, 4)$. Further, let $y\ ϵ\ (7, 8, 9)$. So some of the possible inputs to $f(x,y)$ could be $f(1, 7),\ f(1, 8)$ and so on. Suppose you wanted all possible combinations of x and y in form of a grid. Suppose you represent x, y and (x,y) as shown in Figure 18.4 (Of the book).
The following script shows how the np.meshgrid() method works:
This script is available on page 459 of the book

In [19]:
import numpy as np
X = np.linspace(1, 4, 4) #X-> [1. 2. 3. 4.]
Y = np.linspace(7, 9, 3) #Y-> [7. 8. 9.]
Xgrid, Ygrid = np.meshgrid(X, Y)
print('X->', X,'Y->', Y)
print('Xgrid->', Xgrid)
print('Ygrid->', Ygrid)

X-> [1. 2. 3. 4.] Y-> [7. 8. 9.]
Xgrid-> [[1. 2. 3. 4.]
 [1. 2. 3. 4.]
 [1. 2. 3. 4.]]
Ygrid-> [[7. 7. 7. 7.]
 [8. 8. 8. 8.]
 [9. 9. 9. 9.]]

18.2.9. Using NumPy, SciPy for getting some basic information about a matrix
The SciPy library has a package called linalg for linear algebra. (NumPy also has a linalg package, but it is preferable to use SciPy). NumPy and SciPy have a number of built-in functions to get some basic characteristics of a matrix. For example, you can get characteristics like trace, rank, conjugate, norm, transpose and conjugate transpose. Table 18.2 shows these operations. Table 18.2: Some important methods of the linalg package of NumPy
Functions/ Methods Purpose

1.np.trace(A) Gives trace of A., i.e., sum of all its diagonal elements

  1. linalg.inv(A) Gives inverse. The inverse of a matrix has the property that if it is multiplied by the original matrix, you get an identity matrix
  2. A.T Transpose of matrix
  3. A.H Conjugate transpose
  4. linalg.det(A) Determinant of a matrix
  5. linalg.matrix_norm(A) Gives Frobenius norm (Default)
  6. linalg.norm(A, 1) Gives L1 norm
  7. np.linalg.matrix_rank(A) Rank of matrix. Note function matrix_rank() is available in linalg module of NumPy and not in the linalg module of SciPy

The use of above functions/ methods are shown as follows:
This script is available on page 460 of the book

In [20]:
from scipy import linalg
import numpy as np
np.set_printoptions(precision=2, suppress = True)
A = np.array([[0, -3, -2], 
           [1, -4, -2],
           [-3, 4, 1]])
B = A.T
print('Transpose of A->', B)
C = linalg.inv(A)
print('Inverse of A->', C)
D = linalg.det(A)
print('Determinant of A->', D)
E = linalg.norm(A)
print('Frobenius norm of A->', E)
F = linalg.norm(A)
print('L1 norm of A->', F)
# Note that linalg is available both in numpy and scipy
# Method matrix_rank() is available only in numpy
G = np.linalg.matrix_rank(A)
print('Rank of A->', G)

Transpose of A-> [[ 0  1 -3]
 [-3 -4  4]
 [-2 -2  1]]
Inverse of A-> [[ 4. -5. -2.]
 [ 5. -6. -2.]
 [-8.  9.  3.]]
Determinant of A-> 0.9999999999999991
Frobenius norm of A-> 7.745966692414834
L1 norm of A-> 7.745966692414834
Rank of A-> 3

NumPy provides a function matmul to find the product of two matrices. You know that the product of a matrix with its inverse is an identity matrix. Matrix B is inverse of a square matrix A if $AB = I$ where $I$ is the identity matrix. In common notation, you can also write $B = A^{-1}$. The following code first creates the inverse of a matrix and then does matrix multiplication of the matrix with its inverse:
This script is available on page 461 of the book

In [21]:
import numpy as np
from scipy import linalg
np.set_printoptions(precision=2, suppress = True)
A = np.array([
            [7, 8, 2],
            [3, 5, 1],
            [6, 10, 8]])
B = linalg.inv(A)
print(B)
#np.matmul(B,A) is for multiplying 2 matrices
C = np.matmul(B, A)
# Cast floats to int
print(C.astype(int))

[[ 0.45 -0.67 -0.03]
 [-0.27  0.67 -0.02]
 [ 0.   -0.33  0.17]]
[[0 0 0]
 [0 1 0]
 [0 0 1]]

18.3.1. Solving linear system of equations

You can solve the following system of linear equations using the linalg module of NumPy.
$2x + 3y + z = 13$
$x – y + 2z = 7$
$3x + 4y + z = 22$
The following script shows how to solve these equations:-
This script is available on page 462 of the book

In [22]:
from scipy import linalg
import numpy as np
np.set_printoptions(precision=2, suppress = True)
#2x + 3y + z = 13, x –  y + 2z = 7, 3x + 4y + z = 22
A = np.array([[2, 3, 1], 
           [1, -1, 2],
           [3, 4, 1]])
b = np.array([13, 7, 22])
# Use np.linalg.solve(A, b)
# Return is an array of same shape as b
# Here b is a 1-D array so return will be a 1-D array
x = np.linalg.solve(A, b)
print('x->', x)
# Check solution ie A.x
B = np.dot(A, x)
print('b->', B)
print('Is A.x == b?', np.allclose(B, b))

x-> [11. -2. -3.]
b-> [13.  7. 22.]
Is A.x == b? True

18.3.7. Eigen decomposition
This topic is given on page 465 of the book
The example in the book gives the following matrix and does its eigen decomposition:-

$M = \begin{bmatrix} 3 & 1 & -1 \\[0.3em] 1 & 3 & -1 \\[0.3em] -1 & -1 & 5 \end{bmatrix}$
The script to implement the eigen decomposition is:-
This script is available on page 467-468 of the book

In [23]:
import numpy as np
from numpy import dot #  Otherwise you have to use np.dot
from numpy.linalg import eig, inv
np.set_printoptions(precision=2, suppress = True)
#1--------define a matrix
A = np.array([[3, 1, -1], 
           [1, 3, -1],
           [-1, -1, 5]])
print('A->')
print(A)# Check A
#2(a)---get eigen values as a 1-D array and eigen vectors as matrix Q
eig_arr, Q = eig(A)
print('Eigen values->', eig_arr)# eigen values are a 1-D array
print('Eigen vectors as a matrix->')
print(Q)# eigen vectors are columns of a square matrix
#2(b)------ Convert 1-D eigenvalue array eig_arr into matrix L
L = np.diag(eig_arr)
print('Eigen values as a diagonal matrix->')
print(L)# Confirm L is a diagonal matrix
#3(a)----------Get inverse of the eigen vector matrix Q
Q_inv = inv(Q)
print('Inverse of eigen vector matrix->')
print(Q_inv)# Confirm that inverse of eigen vector matrix is OK
#3(b)----Intermediate result ie B = L.Q_inv
B = L.dot(Q_inv)
print('Intermediate result B->')
print(B) # B is just an intermediate result
#3(c)------ get C = (Q).((L).(Q_inv)) ie C = Q.B
C = Q.dot(B)
print('C = (Q).((L).(Q_inv))->')
print(C)
#4----------Check if A and C are same/ similar
print('Is A similar to C?', np.allclose(A, C))

A->
[[ 3  1 -1]
 [ 1  3 -1]
 [-1 -1  5]]
Eigen values-> [6. 2. 3.]
Eigen vectors as a matrix->
[[-0.41 -0.71  0.58]
 [-0.41  0.71  0.58]
 [ 0.82 -0.    0.58]]
Eigen values as a diagonal matrix->
[[6. 0. 0.]
 [0. 2. 0.]
 [0. 0. 3.]]
Inverse of eigen vector matrix->
[[-0.41 -0.41  0.82]
 [-0.71  0.71  0.  ]
 [ 0.58  0.58  0.58]]
Intermediate result B->
[[-2.45 -2.45  4.9 ]
 [-1.41  1.41  0.  ]
 [ 1.73  1.73  1.73]]
C = (Q).((L).(Q_inv))->
[[ 3.  1. -1.]
 [ 1.  3. -1.]
 [-1. -1.  5.]]
Is A similar to C? True

18.3.8. Singular Value Decomposition
This topic is given on page 469 of the book
The book uses the following $3 \times 2$ matrix and does its SVD

$M = \begin{bmatrix} 1 & -1 \\[0.3em] 0 & 1 \\[0.3em] 1 & 0 \end{bmatrix}$

The script to implement this is as follows:-
This script is available on page 470 of the book

In [24]:
from scipy import linalg
import numpy as np
np.set_printoptions(precision=2, suppress = True)
A = np.array([[1, -1], # A is m x n where m = 3, n = 2
              [0, 1],
              [1, 0]])
# SVD of A ie A = U.S.Vt
U, S, Vt = linalg.svd(A)
# print matrices U, S and Vt and also their shapes
print('U shape->', U.shape, 'U->', U)
print('S shape->', S.shape, 'S->', S)
print('Vt shape->', Vt.shape, 'Vt->', Vt)

U shape-> (3, 3) U-> [[-0.82  0.   -0.58]
 [ 0.41 -0.71 -0.58]
 [-0.41 -0.71  0.58]]
S shape-> (2,) S-> [1.73 1.  ]
Vt shape-> (2, 2) Vt-> [[-0.71  0.71]
 [-0.71 -0.71]]

The remarkable thing to note is that even though S is actually supposed to be a $3 \times 2$ matrix, its shape is shown as $(2,)$ meaning it is a $1D$ array.
This is because for $S$, the only elements which matter are the diagonals, since the rest are all 0. So you can represent the matrix $S$ as a 1D array.
Now suppose you wanted to reconstruct the original matrix A from these three matrices $U$, $S$ and $Vt$. Then you will have to construct a $3 \times 2$ matrix from $S$ first. The number of columns in $U$ are same as number of rows in $S$ and the number of columns in $S$ are same as the number of rows in $Vt$.
Further note that you can get the number of rows in $A$ by using A.shape[0] and number of columns in A by A.shape[1].
You may reconstruct the matrix S from the 1D array S and then do multiplication $U \times S \times Vt$ and get back $A$.
The code is as follows:
This script is available on page 471 of the book

In [25]:
from scipy import linalg
import numpy as np
from numpy import array, dot, shape, diag, zeros
A = array([[1, -1], # A is m x n where m = 3, n = 2
              [0, 1],
              [1, 0]])
# SVD of A ie A = U.S.Vt
U, S, Vt = linalg.svd(A)
# Convert 1-D array S into a 3 x 2 matrix Sm
M,N = A.shape
Sm = linalg.diagsvd(S,M,N)
print(Sm)# Confirm that S converted to 3 x 2 matrix
# Matrix B is product of Sm and Vt. B is an intermediate result
B = Sm.dot(Vt)
print('B = Sm.Vt->', B)
# C is product of U, Sm and Vt
C = U.dot(B)
print('C->', C)
# Check if A and C are same/ similar
print('Is A similar to C?', np.allclose(A, C))

[[1.73 0.  ]
 [0.   1.  ]
 [0.   0.  ]]
B = Sm.Vt-> [[-1.22  1.22]
 [-0.71 -0.71]
 [ 0.    0.  ]]
C-> [[ 1. -1.]
 [-0.  1.]
 [ 1.  0.]]
Is A similar to C? True

18.3.9. Using SVD for dimension reduction
This topic is given on pages 471-473 of the book
In the book dimension reduction is explained with help of an example of 4 readers categorizing 6 books.
This creates a $4 \times 6$ matrix as shown below:-
$M = \begin{bmatrix} 10 & 8 & 9 & 1 & 0 & 0 \\[0.3em] 9 & 10 & 9 & 0 & 1 & 1 \\[0.3em] 1 & 0 & 1 & 10 & 8 & 1 \\[0.3em] 0 & 1 & 1 & 9 & 8 & 0 \end{bmatrix}$
After doing the categorization, you get:-
$Ʃ = [22.81, \ 17.26, \ 2.14, \ 0.96]. $
Out of these 4 values, only the first and second are significant. So you can think of these two as the two categories in which the books are falling into. The 3rd and 4th entries can be thought of as noise. Keep only the first 2 entries of Ʃ and convert the other 2 to 0.
Call this modified $Ʃ$ as $Ʃmod$. Then we have:-
$Ʃmod = \begin{bmatrix} 22.81 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & 17.26 & 0 & 0 & 0 & 0 \end{bmatrix}$

So only 2 rows and 2 columns are kept. Then we have:-
$Ʃmod = \begin{bmatrix} 22.81 & 0 \\[0.3em] 0 & 17.26 \end{bmatrix}$

The following code repeats all the steps outlined in the book (It also uses a module linalg for linear algebra from the SciPy library):
This script is available on page 474-475 of the book

In [26]:
from scipy import linalg
from numpy import array, dot, shape, diag, zeros
A = array([[10, 8, 9, 1, 0, 0],
           [9, 10, 9, 0, 1, 1],
           [1, 0, 1, 10, 8, 1],
           [0, 1, 1, 9, 8, 0],])
# 1---Create SVD of matrix A and convert S from 1_d array to matrix---
#1(a)---SVD of A ie A = U.S.Vt
U, S, Vt = linalg.svd(A)
print('U->', U)
print('S->', S)
print('Vt->', Vt)
print(type(S))
# 1(b) ---Create Sm ie matrix from linear 1-D array S----
M,N = A.shape
Sm = linalg.diagsvd(S,M,N)
print('Sm->', Sm)
#2 Modify Vt to have 2 rows. Modify S to be 2x2. Modify Y to have 2 columns
S2 = Sm[:2, :2]
Vt2 = Vt[:2, :]
U2 = U[:, :2]
print('Vt2->', Vt2) # Confirm Vt2 has 2 rows
print('S2->', S2) # Confirm S2 is 2 x 2
print('U2->', U2)# Confirm U2 has 2 columns
#3------- Reconstruct Y = U2.S2.Vt2 ---
#3(a)----X is an intermediate dot product S2.Vt2
X = S2.dot(Vt2)
#3(b)-----Y = ((U2.(S2.Vt2)) ie Y = U2.X -----
Y = U2.dot(X)
print('Y->', Y)
#4---Compare original matrix A to matrix Y using Frobenius norm
Frob_Y = linalg.norm(Y)
print('Frobenius norm of Y->', Frob_Y)
Frob_A = linalg.norm(A)
print('Frobenius norm of A->', Frob_A)

U-> [[-0.67  0.17  0.64  0.33]
 [-0.7   0.19 -0.62 -0.32]
 [-0.18 -0.71  0.31 -0.61]
 [-0.17 -0.66 -0.33  0.65]]
S-> [22.81 17.26  2.14  0.96]
Vt-> [[-0.58 -0.55 -0.56 -0.18 -0.16 -0.04]
 [ 0.16  0.15  0.11 -0.74 -0.62 -0.03]
 [ 0.54 -0.65  0.09  0.33 -0.39 -0.15]
 [-0.19  0.11  0.14  0.05 -0.02 -0.96]
 [ 0.04 -0.49  0.45 -0.49  0.56 -0.03]
 [-0.56 -0.09  0.67  0.25 -0.35  0.22]]
<class 'numpy.ndarray'>
Sm-> [[22.81  0.    0.    0.    0.    0.  ]
 [ 0.   17.26  0.    0.    0.    0.  ]
 [ 0.    0.    2.14  0.    0.    0.  ]
 [ 0.    0.    0.    0.96  0.    0.  ]]
Vt2-> [[-0.58 -0.55 -0.56 -0.18 -0.16 -0.04]
 [ 0.16  0.15  0.11 -0.74 -0.62 -0.03]]
S2-> [[22.81  0.  ]
 [ 0.   17.26]]
U2-> [[-0.67  0.17]
 [-0.7   0.19]
 [-0.18 -0.71]
 [-0.17 -0.66]]
Y-> [[9.31 8.85 8.84 0.53 0.54 0.5 ]
 [9.66 9.18 9.16 0.45 0.48 0.52]
 [0.53 0.49 1.03 9.82 8.25 0.53]
 [0.51 0.47 0.97 9.2  7.73 0.5 ]]
Frobenius norm of Y-> 28.609288664939417
Frobenius norm of A-> 28.705400188814647

Assignment
Using SVD for image compression
This topic is given on page 478-479 of the book
(This section uses the Matplotlib library. If you are not familiar with it, first read the chapter on Matplotlib.)
In the chapter, the SVD algorithm was discussed in details in context of Dimension reduction. This assignment is about actually using the SVD algorithm for dimension reduction of an image leading to image compression. For dimension reduction of an image, you need to first convert it to a matrix. Then you do dimension reduction on the matrix and then convert the matrix back to an image. The resultant image will be compressed. For doing this, the Image module of the PIL/ Pillow library is used. The steps are:

  • Load an image from your hard disk using the Image.open(path_name) method.
  • Get size of the image in terms of its width and height.
  • Create a Matplotlib figure. On this figure, create three sub-plots. On first sub-plot you will generate the colour uncompressed image. On second sub-plot, you will generate uncompressed grey (black and white) image. On the third sub-plot you will generate compressed grey image.
  • However, before you can do dimension reduction, you need to convert the image into a matrix.
  • The steps are:
    (a) convert image into sequence using the getdata() method of a PIL image,
    (b) convert the sequence into a Python list,
    (c) convert the list into a numpy.ndarray,
    (d) reshape the numpy.ndarray into same 2D as original image whose width and height had been saved into two variables named width and height and
    (e) finally convert the reshaped ndarray into a NumPy matrix.
  • Apply SVD on the resultant matrix like any other matrix. In the present case, you have reduced S to 1/20th of original. Accordingly also reduce the columns of U and the rows of Vt.
    The script is as follows:
In [27]:
% matplotlib inline
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
path = r'H:\=Python\test-image\view2.jpg'# Path to your file here
#1--- Open image in PIL. Get width and height of image
img = Image.open(path)
#details of size parameter of PIL image at link in next line
#https://pillow.readthedocs.io/en/4.1.x/reference/Image.html#PIL.Image.size
width, height = img.size
print('Image width->', width, 'Image height->', height)#width, height of image
#2---- Create matplotlib figure for plotting.
fig = plt.figure(figsize=(8, 10))
#3-------- Plot uncompressed color image
ax_col = fig.add_subplot(131)
ax_col.imshow(img)
#4------ Plot uncompressed grey image
ax_grey = fig.add_subplot(132)
imggray = img.convert('LA')
ax_grey.imshow(imggray, cmap = 'gray')
#5---Convert image data into a 2-D matrix
# For getdata() see:-
# http://effbot.org/imagingbook/image.htm#tag-Image.Image.getdata
# im.getdata() Returns the contents of an image as a sequence
img_seq = imggray.getdata(band = 0)
#To convert it to an ordinary sequence use list(im.getdata()).
img_list = list(img_seq)
img_arr = np.array( img_list, float)
print('shape array->', img_arr.shape)
# The shape property can be used to get the shape as well as to set the shape.
# Here used for setting the new shape.See following link
#https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.shape.html
img_arr.shape = (height, width)
print('shape array after ', img_arr.shape)
#Convert ndarray to matrix
img_mat = np.matrix(img_arr)
#6----SVD and dimension reduction ie compression of image
U, S, Vt = np.linalg.svd(img_mat) #single value decomposition
# Reduce S to 1/20th of original
c1 = int(len(S)/20)
print('c1->', c1)
# Reduce columns of U, make diagonal matrix S and reduce rows of Vt
img_c1 = np.matrix(U[:, :c1]) * np.diag(S[:c1]) * np.matrix(Vt[:c1,:])
ax_c1 = fig.add_subplot(133)
plt.imshow(img_c1, cmap = 'gray')

Image width-> 862 Image height-> 511
shape array-> (440482,)
shape array after  (511, 862)
c1-> 25
Out[27]:
<matplotlib.image.AxesImage at 0x8637410>